#!/usr/bin/python
import sys
import os

# Read data.
def read_file(filepath):
	# define data structure.
	data = {}
	
	# Check if zipped.
	zipped = False
	if filepath.count(".gz"):
		zipped = True
	
	# Open file.
	seq = ""
	seqheader = ""
	first = True
	if zipped == True:
		fin = gzip.open(filepath,"r")
	else:
		fin = open(filepath,"r")
		
	# Read the whole thing for speed.
	lines = fin.readlines()
	fin.close()
	
	# Loop over read lines.
	for line in lines:
	#cnt = 0
	#for line in fin:
		# Skip comments.	
		if line[0] == "#": continue
	#	cnt += 1
	#	if cnt > 100000: break


		# first time.
		if first == True and line[0] == ">":
			seq = ""
			seqheader = line.strip().replace(">","")
			first = False
			continue		

		# finished seq.
		if line[0] == ">":
			# save data.
			data[seqheader] = seq.strip()

			# Reset variables.
			seq = ""
			seqheader = line.strip().replace(">","")
			continue

		# get seq.
		seq += line.strip()
			
	# Handle last close and return.
	data[seqheader] = seq.strip()
	return data

# Read ref.
print "reading seqs"
ref = read_file(sys.argv[2])
print "done"

# read maps.
print "reading maps"
fin = open(sys.argv[1],"rb")
lines = fin.readlines()
fin.close()
print "done"
print "we have to map,", str(len(lines))

# open output.
fout = open(sys.argv[3],"wb")

# Check mapping
cnt = 0
for line in lines:
	cnt += 1
	if cnt % 100000 == 0:
		print cnt
	
	# split
	tmp = line.split("\t")
	
	# get matchs.
	i = tmp[5].find("M")
	
	# get first hard clip.
	j = tmp[5].find("H")
	
	# Get match.
	try:
		match = int( tmp[5][j+1:i] )
	except:
		continue
	
	# Get percentage alignment.
	rname = tmp[0].replace(">","")
	if rname in ref:
		
		# Apply filter.
		id = float(match) / float(len(ref[rname]))
		if id >= .9:
			fout.write( line )
			
fout.close()
	

